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We investigate a gas of superfluid fermionic atoms trapped in two hyperfine states by a spherical 
harmonic potential. We propose a new regularization method to remove the ultraviolet divergence 
in the Hartree-Fock-Bogoliubov equations caused by the use of a zero-range atom-atom interaction. 
Compared with a method used in the literature, our method is simpler and has improved convergence 
properties. Then we compare Hartree-Fock-Bogoliubov calculations with the semiclassical local- 
density approximation. We observe that for systems containing a small number of atoms shell effects, 
which cannot be reproduced by the semiclassical calculation, are very important. For systems with 
a large number of atoms at zero temperature the two calculations are in quite good agreement, 
which, however, is deteriorated at non-zero temperature, especially near the critical temperature. 
In this case the different behavior can be explained within the Ginzburg-Landau theory. 
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I. INTRODUCTION 

In the last few years an increasing interest has 
been directed towards ultracold gases of trapped 
fermionic atoms. Many experimental efforts are 
made to develop and improve the techniques for 
trapping and cooling fermionic atoms like, for in- 
stance, 40 K and 6 Li. An interesting aspect of 
trapped fermionic atoms in comparison with other 
Fermi systems is that parameters such as the tem- 
perature, the density, the number of particles, and 
even the interaction strength are tunable exper- 
imentally. By tuning the magnetic field in the 
vicinity of a Feshbach resonance , the scattering 
length, which is related to the interaction strength, 
can be changed. This offers a wide range of pos- 
sibilities to investigate the behaviour of these sys- 
tems in different experimental conditions. By us- 
ing optical or magnetic traps, temperatures of 
about jTp have been achieved 0, 0, 0], where 
Tp = ep/ks is the Fermi temperature. 

All these efforts are mainly directed to the real- 
ization and detection of a phase transition to the 
superfluid phase below some critical temperature 
Tc- In order to have a s-wave attractive inter- 
action among the atoms, which can give rise to 
s-wave pairing correlations below Tc, the atoms 
have to be trapped and cooled in two different hy- 
perfine states. This has been achieved in a recent 
experiment 5] , where also the Feshbach resonance 
in the 6 Li scattering amplitude has been used to 
enhance the scattering length. It seems that in the 
same experiment some signals indicating a super- 
fluid phase transition have been observed. 

From the theoretical point of view many calcula- 
tions have been performed in order to predict and 
study the equilibrium properties of the trapped 
system when the phase transition takes place. So 
far all these calculations are based on the mean- 



field approach. In Ref. Q the trapped Fermi gas 
was treated in local-density approximation (LDA) , 
where the system is locally treated as infinite and 
homogeneous. In Ref. some corrections to the 
LDA for temperatures near Tq were obtained in 
the framework of the Ginzburg-Landau (GL) the- 
ory. The first approach fully taking into account 
the finite system size was introduced in Ref. 
and studied further in Refs. 0, EJ- It consists 
in a Hartree-Fock-Bogoliubov (HFB) calculation, 
analogous to calculations done in nuclear physics, 
where the mean field and the pairing properties of 
the system are treated self-consistently. In Ref. Q 
also a regularization prescription for the pairing 
field was developed: Since the densities in the traps 
are very low, the atom-atom interaction can be ap- 
proximated by a zero-range interaction. However, 
this leads to an unphysical ultraviolet divergence 
of pairing correlations which has to be removed. 

In spite of the possibility to perform full HFB 
calculations, it should be mentioned that these cal- 
culations are numerically very heavy and therefore 
limited to moderate numbers of particles. Another 
shortcoming of present HFB calculations is that 
they are restricted to the case of spherical sym- 
metry, while the traps used in the experiments 
are usually strongly deformed. Hence, to describe 
trapped systems under realistic conditions, one has 
to rely on calculations within the LDA. This is a 
quite embarrassing situation, since even for large 
numbers of particles the results of HFB and LDA 
calculations have not always been in good agree- 
ment (see results shown in Ref. y|). 

In this paper we will present a detailed com- 
parison between HFB and LDA calculations. In 
particular, we will show that the disagreement be- 
tween HFB and LDA calculations which has been 
found in Ref. is to a certain extent caused by 
the use of an unsuitable regularization prescription 
for the pairing field in the HFB calculations. We 
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will present a modified regularization prescription 
which was originally developed for HFB calcula- 
tions in nuclear physics and which is much 
easier to implement numerically. (As we learned 
after sending the first version of our manuscript, 
Nygaard et al. used the same prescription in their 
calculation of a vortex line in a dilute superfluid 
Fermi gas 12], without giving a description of this 
scheme.) Due to its improved convergence proper- 
ties, this scheme leads to more precise results for 
the pairing field, which in the case of large num- 
bers of atoms agree rather well with the results 
of the LDA at least at zero temperature. At non- 
zero temperature, however, the differences between 
HFB and LDA results turn out to be important 
even for very large numbers of particles. For ex- 
ample, we find that the critical temperature ob- 
tained within the LDA is too high, and that the 
pairing field profile near the critical temperature 
is not well described by a LDA calculation: we 
show with the HFB approach that it actually has 
a Gaussian shape, as it was predicted in the frame- 
work of the GL theory in Ref. 0. 

The paper is organized as follows: In Sec. ILTlwe 
will present the adopted formalism with a partic- 
ular attention on the description of the regular- 
ization techniques. In Sec. IIIII we will show some 
comparisons between HFB and LDA calculations 
and illustrations of the results obtained with dif- 
ferent choices for the regularization method. We 
will also discuss results obtained for non-zero tem- 
peratures and verify the quantitative predictions 
of the GL theory. Finally, in Sec. |IV| we will draw 
our conclusions. 



II. THE FORMALISM 

In this paper we will consider a spherically sym- 
metric harmonic trap with trapping frequency w, 
where N atoms of mass m populate equally two 
different spin states f and J., i.e., Nf = N±. As 
mentioned in the introduction, the low density of 
the system allows to introduce a contact interac- 
tion for the atoms, caracterized by the s-wave scat- 
tering length a. The hamiltonian reads 

H = T + ± *™^ + ^ E^-r,-), (1) 

3=1 i<3 

where T is the kinetic term. For convenience let 
us introduce a coupling constant g defined as: 



Since we are considering attractive interactions, we 
have a < and, consequently, g < 0. To simplify 



the notation, we will use in what follows the "trap 
units" , i.e. 

m = u> = h = ks = 1 • (3) 

Thus, energies will be measured in units of hui, 
le ngths in units of the oscillator length l\ lo = 
yj h/ (mu)), and temperatures in units of hj/ks- 

Before describing the HFB approach, let us add 
some comments on the validity of the hamilto- 
nian The parametrization of the interaction 
in terms of the free-space s-wave scattering length 
a is valid at very low densities, where the distance 
between particles is much larger than \a\. However, 
if the distance between particles becomes compara- 
ble with \a\, the bare interaction has to be replaced 
by a density-dependent effective interaction, as it 
is done in nuclear physics (see also 0). This is 
particularly important in the vicinity of a Fcsh- 
bach resonance, where \a\ becomes very large. In 
this case it might be necessary to include the Fes- 
hbach resonance as a new degree of freedom into 
the Hamiltonian [bH ]. 

A. HFB approach and regularization 
procedure 

The hamiltonian (|TJ will be treated within the 
mean-field approximation. We will not go into de- 
tails here as the formalism has been introduced and 
extensively illustrated in Ref. The Hartree- 
Fock-Bogoliubov (HFB) or Bogoliubov-de Gennes 
[TH Il6l ] equations read: 

[H + W(R)]u a (R) + A{R)v a (R) = E a u a {R) , 
A(R)u a (R) - [H + W(R)]v a (R) = E a v a (R) , 

(4) 

where a collects all quantum numbers except spin 
(n, I, m), u a and v a are the two components of the 
quasiparticle wavefunction associated to the en- 
ergy E a , and Hq is the following single-particle 
hamiltonian: 

H = T + U -ti, (5) 

where Uq = ^r 2 is the harmonic trapping potential 
and n the chemical potential. The Hartree-Fock 
mean field W(R) in Eq. is expressed by 

W(R)=g^2{[v a (R)[ 2 [l-f(E a )] 

(X 

+ \u a (R)\ 2 f(E a )} , (6) 
where f{E a ) is the Fermi function: 

f( E «) - ^JT-^ ■ ( 7 ) 
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With a zero-range interaction the pairing field 
A(i?) appearing in Eq. (0} would usually be de- 
fined as A(R) = -g(^>i(R)^i(R)), where is 
the field operator associated with the spin states J. 
and f. However, this expression is divergent and 
must be regularized. The regularization prescrip- 
tion proposed in Ref. @ consists in using the pseu- 
dopotential prescription |l7j : 

A(R) = -gVm^r^R+^^R- §))] . 

(8) 

In practice, Eq. (JSJ is evaluated as follows: It 
is possible to show that the expectation value 
(*t(-R + r/2)*j.(fl - r/2)) diverges as A/(4nr) 
when r — > if a zero-range interaction is used. 



Now one adds and subtracts from this expecta- 
tion value the quantity ^A(i?)G° (H, r), where G° 
is the Green's function associated to the single- 
particle hamiltonian Ho, Eq. JSJ, and calculated 
for the chemical potential \x: 

where 0° denotes the eigenfunction of f/o with 
eigenvalue e° — /i. One can demonstrate that this 
Green's function diverges as l/(27rr) when r — ► 0. 
Expressing (Wf^j,) in terms of the wave functions 
u and v, one can write the pairing field A as 



A(R) 



a 

A(i?) 0°(Jl+|)^(fl-? 



(10) 



The sum over a is no longer divergent for r — ► 0, 
since the divergent part of — lAG° cancels the di- 
vergent part of (vE'f^I'x). Thus, we can take the 
limit r — > of this sum. On the other hand, the 
divergence of the last term is removed by the pseu- 
dopotential prescription, which selects only the 
regular part of the Green's function G°: 



l™l r [rGl(R>r)]=G^(R). 
Finally, A can be expressed as follows: 



(11) 



A(R) =-gJ2 (u a (RK(R) [1 - 2f(E a )} 

(12) 



A* 

Once the regular part of the Green's function is 
calculated for a given chemical potential /u, |8( , the 
HFB equations are solved self-consistently. 

In practice, it is of course impossible to extend 
the sum over all states a and one has to intro- 
duce some cutoff. However, since the sum over a 
converges, the cutoff should not affect the results 
if it is chosen sufficiently high. We will discuss 
about the rapidity of convergence of the regular- 
ization procedure presented here with respect to 
the introduced energy cutoff. We will show that 
the convergence is quite slow. Moreover, the cal- 
culations can become heavy when systems with a 



large number of atoms are treated, as the function 
QQreg ^ j-^ ca j ciua ted for a large value of the 
chemical potential. A way to simplify the regular- 
ization procedure and to avoid to calculate G° res 
is proposed in Ref. ^lj , where the procedure of Q 
is extended to calculations for nuclear systems. We 
will describe this method in next subsection. 



B. Thomas-Fermi approximation in the 
regularization procedure 

In Ref. p"H a simpler regularization procedure 
was proposed where the Thomas-Fermi approxi- 
mation (TEA) is used to calculate the regular part 
of the Green's function. To that end let us write 
the Green's function G° by adopting the TFA for 
the sum over the states corresponding to oscilla- 
tor energies above some suffiently large value 



G°(fl,r)« J2 



nlm 



d 3 k 



, (13) 



where 



k c (R) = V^ec - R 2 = \/2N c + 3-R 2 . (14) 
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Observing that and using Eq. (|13fl . we can write the regular part 

of the Green's function as follows: 

" + °° d 3 k e tk r _ 1 

J^Kf if ~ 2^ ( ' 



nhn 



^ f _ ,. + / / (9.^3 I fc 2 i -R 2 ,, fc 2 j / Ott)3 k 2 ■ \ ®> 

nlm € nl M Jk c (R){^) ^ — + — - ^ — 7 ^0 W — 

Evaluating the integrals over k and summing over the magnetic quantum number m, we obtain 

pO^M- V ( 2 * + 1 Kt(r) | fcg-W, fcc(r) + fcg.(r) fc c (r) 
^ 2^ 47r(e°;-/x) 2^ 2 fe c (r) - fe° (r) - 2 ' 



nl 
e1,<e c 



where R n i are the radial parts of the oscillator wave 
functions and 



k° F (r) = y/2fi - r 2 



(18) 



is the local Fermi momentum. As noted in 
Ref . ' this method can be used beyond the clas- 
sical turning point (characterized by k F (r) = 0) 
by allowing for imaginary values of k F (r). The 
case that kc(r) becomes imaginary will not be 
considered, because we assume that Nc is suffi- 
ciently large such that the pairing field can be ne- 
glected in the regions where kc{r) is imaginary. It 
should also be pointed out that already for, say, 
Nc > jtx + 10, Eq. 1)17(1 is an extremely accurate 
approximation to G" res , and gives results which 



J 



are almost undistinguishible from those obtained 
by the numerically heavy algorithm proposed in 
Ref. @. 

Now let us substitute Eq. JT7|) into Eq. JT^J). We 

have to choose a cutoff for the sum over single- 
particle states. Instead of choosing a cutoff for 
the quasiparticle energies E a , as it is done in Ref. 
[ll| . we can likewise restrict the sum in Eq. (|12|) 
to the states corresponding to those appearing in 
the sum in Eq. (|17|l . This is the natural choice if 
one obtains the wave- functions u a and v a and the 
quasiparticle energies E a by solving Eq. (0J in a 
truncated harmonic oscillator basis containing the 
states satisfying e° ; < ec = Nc + §• In this way 
we obtain the following simple formula for the gap: 



»,s 2Z + 1 , . . , r „ , n A(r) /fc£(r) , k c (r) + k° F (r) k c (r)\ 

Mr) = ~9 J2 -^(rK.(r) [I - 2f(E nl) ] ~ 9^{^ ^ _ ^ - ■ (19) 

Finally, this can be rewritten in terms of a position dependent effective coupling constant: 



^{r)=-g eff (r) J2 



21 + 1 

47T 



u n i{r)v nl {r) [1 - 2f(E nl )} , 



(20) 



where 



1 , 1 (k° F {r). k c {r) + k° F {r) 



1 1 1 / 

ffir) q 2tt 2 \ 



g eff (r) g ' 2ir 2 K 2 hl k c ( r ) - k F (r) M0 



(21) 



We stress again that the results obtained with prescription (a) , coincide with the results obtained 
this regularization prescription, from now on called 
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with the prescription introduced in Rcf. 8] . 

However, it will turn out that it is useful to in- 
troduce the following modification of the method: 
Let us replace everywhere k F (r) by the local Fermi 
momentum taking into account the full potential 
(trapping potential Uq plus Hartree-Fock potential 
W): 

k F (r) = ^2(1 - r 2 - 2W(r) . (22) 

Formally this replacement does not change any- 
thing: Instead of adding and subtracting the term 
|A(i?)G° (R, r) from the divergent expectation 
value (* T (J2+r/2)* i ( J R-r/2)) with G° being the 
Green's function corresponding to the harmonic 
oscillator potential Uq, we can also add and sub- 
tract a similar term involving the Green's function 
G M corresponding to the full potential Uq + W. 
Also from Eq. (|21|l it is evident that in the limit 
Nc — > oo [i.e., kc{r) — * oo] the results will be inde- 
pendent of the choice of kp . However, we will see 
that the convergence of this modified scheme, from 
now on referred to as scheme (b) , is very much im- 
proved. Thus, it is possible to use a much smaller 
cutoff Nc without having a strong cutoff depen- 
dence of the results. 



C. Local-density approximation 



If the number of particles becomes very large, it 
is natural to assume that the system can be treated 
locally as infinite matter with a local chemical po- 
tential given by /i — Uq(t). This assumption leads 
directly to the local-density approximation (LDA) . 
Formally, the LDA corresponds to the leading or- 
der of the Wigner-Kirkwood h expansion, which 
is at the same time an expansion in the gradients 
of the potential ^5|. Thus it is the generaliza- 
tion of the standard Thomas-Fermi approximation 
(TFA), which also corresponds to the leading or- 
der of an h or gradient expansion, to the superfluid 
phase. Here we will adopt the name LDA in order 
to avoid confusion with the full HFB calculations 
using the TFA only in the regularization prescrip- 
tion, as discussed in Sec. Ill Til But in the literature 
also the name TFA is adopted. 

In the case of a zero-range interaction, the LDA 
(or TFA) amounts to solving at each point r the 
following non-linear equations for the mean field 
W{r) and the pairing field A(r): 



W(r) = y -p{r) = g 
d 3 k 



d 3 k 
(2^ 



[1 



A(r) = - 5 J -1* ([i- 2 /(25(r,fe))] 



2f(E(r,k))] 
A(r) 



e{r,k) - /i 
2E(r,k) 
A(r) 



2E(r,k) 2(e(r,fc)-/i) 



(23) 
(24) 



where 

e{r,k) = — + U (r) + W(r), (25) 
E(r, k) = v/(e(r,fc)- M ) 2 + A 2 (r) . (26) 

The last term in Eq. 1241) has been introduced in 
order to regularize the ultraviolet divergence. In 
fact, the pseudopotential prescription used in the 
previous subsections was originally motivated by 
the fact that it reduces to such a term if it is ap- 
plied to a homogeneous system 0, ^| . A more rig- 
orous justification of this term is that it appears if 
one renormahzes the scattering amplitude of two 
particles in free space |lgj . 

Let us first consider the case of zero tempera- 
ture, T = 0. In this case, and if the gap A is small 
compared with the local Fermi energy tp — \k 2 F , 
Eqs. l|2"3"|) and can be solved (almost) analyt- 
ically. Under these conditions the density practi- 
cally coincides with the density obtained for A = 0, 
where Eqs. (|23|) . (|25|1 . and H2(jfl can be transformed 



into a cubic equation for the local Fermi momen- 
tum: 

g kM + k^)_ + Uo{r) _ fi = (27) 

For a given local Fermi momentum and under 
the assumption that corrections of higher order in 
A/ep are negligible, Eq. (|24|) can be solved ana- 
lytically. The result is the well-known formula 

A(r ) = 8e,(r-)exp(-2-^L_). (28) 

Now we turn to the case of non-zero tempera- 
ture, but we want to consider only temperatures 
below the critical temperature, i.e., < T < Tc- 
Therefore, we can neglect the influence of the tem- 
perature on the density and have to consider only 
the temperature dependence of A. Let us denote 
the gap at T = by Ao. Then the gap at non- 
zero temperature can be obtained from the approx- 



G 



imate relation H 



-m^W^^KP- (29) 



Ao(r) 



VC 2 + A 2 (r) 



The solution of this equation leads to a universal 
function which gives the ratio A / Ao as a function 
of T/T c , with T c ~ 0.57A . Note that, within the 
LDA, the critical temperature is a local quantity, 
T c = T c {r). 

In order to compare the LDA with the HFB 
theory, with special emphasis on the regulariza- 
tion prescription, wc will now introduce a regular- 
ization scheme for the gap equation within LDA 
which is slightly different from Eq. (|24|l . First of 
all, if we want to investigate the cutoff dependence, 
we have to introduce a cutoff in Eq. (|2"4"j) • Secondly, 
the regularization term introduced in Eq. (12411 



corresponds to the regularization prescription (b) 
described at the end of the previous subsection, 
which is different from that introduced in Ref. || 
and from the regularization scheme (a). If we want 
to compare the LDA results with HFB results ob- 
tained with the original prescription or with the 
prescription (a), which involves the Green's func- 
tion G° of the potential Uq and not the Green's 
function G M of the full potential Uq + W, we have 
to replace the energy e(r, k) appearing in the reg- 
ularization term by 



e°(r,k) = Y + Uo(r) 



(30) 



Thus, the gap equation within LDA suitable for 
comparison with the regularization scheme (a) 
reads 



am = -» jf" ^ (ii - w. m^ggj - 2(t ^ ( ;,-M) ) • (31) 

At zero temperature, T = 0, it is again possible to solve this equation analytically, with the result 



lk c (r)-k F (r) 7T fc°(r) k c (r) - k° F (r) 



A(r) = 8eFir l kc(r) + k F (r) 6XP V 2 2M^)H " W ln kc(r) + k%(r) ) ' (32) 



r 



The result corresponding to the regularization 
scheme (b), Eq. (|28[) . is recovered from this result 
by replacing k F by kp. In this case there is no 
cutoff dependence at all, but one should remem- 
ber that in deriving Eq. (|32[) we have implicitly 
assumed that the cutoff lies above the Fermi sur- 
face. A weak cutoff dependence would appear only 
if corrections to Eq. (|32[) of higher order in A/e F 
were included. 



III. NUMERICAL RESULTS 

In this section we will present some numerical 
results. In particular, we will investigate the con- 
vergence properties of the different renormalization 
methods. Then, we will discuss the validity of the 
LDA at zero temperature. Finally, we will compare 
HFB and LDA calculations at non-zero tempera- 
ture. 

In our numerical calculations we will use for the 
coupling constant the value g = —1 (in units of 
H 2 lho/m). If we consider 6 Li atoms with scatter- 
ing length a — — 2160ao [2(j , where ao = 0.53 A is 
the Bohr radius, this value of g corresponds to a 



trap with oj = 2ir x 817 Hz. (Before relating this 
to real experimental conditions, one should how- 
ever remember that in the experiments the trap is 
usually axially deformed, with a low longitudinal 
trapping frequency lo z and a high transverse trap- 
ping frequency lu±. For example, in the experiment 
described in Ref. [5j , the trapping frequencies were 
given by uj z = 2tt x 230 Hz and ll>j_ = 2ttx 6625 Hz.) 
The choice g = — 1 also facilitates the comparison 
of our results with those from Ref. Q, where the 
same value for g was used. 

A. Convergence of the regularization 
methods 

In this section we will discuss the convergence 
rates with respect to the cutoff used in the numer- 
ical calculations for different choices for the regu- 
larization procedure. As in Sec.[H]we denote by (a) 
the HFB calculations made with the choice of k F 
given by Eq. (|18|) , and and by (b) the calculations 
made with the choice where k F is replaced by k F 
as given by Eq. (|22p. For our comparison we use 
a chemical potential fj, = 32fiio, the corresponding 
number of atoms in the trap is N « 1.7 x 10 4 . 
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In Figs, n and [21 we present the pairing field 
A calculated at zero temperature within the HFB 
and LDA formalisms for different values of the cut- 
off Nc from 50 up to 125. The results shown 
in Fig. ^ have been obtained with the choice (a) 
for the regularization for both the HFB and LDA 
calculations. We verified that the HFB calcula- 
tions with the exact Green's function G° re9 (with- 
out TEA) give practically the same results as the 
method HFB (a) for all the values of the cutoff. 
This means that the TFA in the regularization pro- 
cedure is very satisfying and reproduces well the 
regular part of the oscillator Green's function. 

We observe in Fig.^that the agreement between 
LDA and HFB is reasonable for all values of the 
cutoff Nc- We also notice that for Nc = 125, 
which is the maximum value that we considered, 
the convergence has not yet been reached and 
therefore the pairing field would grow further if 
we could increase the cutoff above 125. In Fig. [21 
we present the same calculations made with the 
choice (b) for the regularization. Remember that 
with this choice, the pairing field within LDA is 
independent of Nc once Nc lies above the Fermi 
surface. On the other hand, the HFB results sat- 
urate quite fast and are already very close to con- 
vergence for Nc — 75. Again, the LDA and HFB 
results are in reasonable agreement. 

By comparing Figs.n an d[21 0ne observes clearly 
that the calculations (a), Fig. ^ are still quite 
far from convergence even for the highest consid- 
ered cutoff. We argue that the convergence rate 
of method (a), which is the same convergence rate 
as that of HFB without TFA in the regulariza- 
tion prescription is much slower than that of 
method (b) . This is more evident in Fig. [3 where 
we plot the HFB values of the pairing field in the 
center of the trap, A(0), for the two regularization 
prescriptions (a) (stars) and (b) (diamonds) as a 
function of the cutoff Nc- We also plot the results 
obtained within the LDA(a) (full line) and LDA(b) 
(dashed line) up to a cutoff of Nc = 10 4 . In the 
inset of the figure we magnify the region of cutoff 
values between 50 and 150. We can observe in the 
inset that the LDA (a) curve fits well the calculated 
points for HFB (a). We noticed that the LDA (a) 
results converge slowly towards a pairing field of 
about 6.86huj, at a very high cutoff, Nc = 10 6 . 
For Nc = 10 3 the pairing field in LDA (a) is still 
only 6.37huj. This very slow convergence rate can 
be understood within the LDA by taking the ratio 
of the pairing fields corresponding to the methods 
(a) and (b). Using Eq. (|32[) in the limit of very 
large kc, one can derive the relation 



^LDA 



\g\V2[»-W(r)} 



, (33) 



a HFB(a) (N c =50) 
x HFB(a) (N c =75) 
+ HFB(a) (N c =100) 
♦ HFB(a) (N c =125) 

— LDA(a) (N c =50) 
LDA(a) (N c =75) 

— LDA(a) (N c =100) 
LDA(a) (N =125) 




FIG. 1: Pairing field A (in units of hu)) as a func- 
tion of the distance r (in units of Iho) from the center 
of the trap, calculated for the parameters = ?>2hu} 
and g = — lft 2 /h /m, corresponding to N ~ 1.7 x 10 4 
particles in the trap. The different curves have been 
obtained within the HFB and LDA formalisms using 
the regularization prescription (a) for different values 
of the cutoff Nc- 



a HFB(b) (N c = 


=50) 


* HFB(b) (N c = 


=75) 


+ HFB(b) (N c = 


=100) 


♦ HFB(b) (N c = 


=125) 


LDA(b) 





< 4 - 



2 - 




A LZM(b) (r) 37r 2 VAc 
where W(r) represents the Hartree field (in the 



FIG. 2: Same as Fig. Q but with regularization pre- 
scription (b). Remember that with this prescription 
the LDA result [Eq. 1321 1 is independent of the cutoff 
N c . 



present case, W(0) ~ — 167uj). 

As the agreement between LDA(a) and HFB (a) 
is good in the region up to Nc — 125, we suppose 
that the convergence rate for HFB (a) is the same as 
for LDA(a). On the contrary, within HFB(b) the 
values of the pairing field in the center of the trap 
are 6Mfuv for N c = 100 and 6.86/iu; for N c = 125: 
we conclude that the convergence in this case is 
much faster. In what follows we will always use 
the method (b) for the regularization procedure. 
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FIG. 3: Value of the pairing field in the center of the 
trap, A(0) (in units of hu>), as a function of the cutoff 
Nc , obtained from HFB calculations with the regular- 
ization methods (a) (stars) and (b) (diamonds), and 
from the LDA, method (a) (solid line) and method (b) 
(dashed line). The parameters fi and g are the same 
as in Fig. 



FIG. 4: Value of the pairing field in the center of the 
trap, A(0) (in units of hu), as a function of the number 
of particles, N, obtained from HFB (solid line) and 
LDA (dashed line) calculations [regularization method 
(b), cutoff No = 100, coupling constant g = — 1 in trap 
units] . 



B. Validity of the LDA at zero temperature 

As mentioned before, the parameters used for 
the calculations shown in Figs. ^ [21 an d El corre- 
spond to a trap with about 1.7 x 10 4 atoms. In this 
case we found a good agreement between the nu- 
merical HFB results and the results obtained from 
the LDA. However, one might wonder under which 
conditions the LDA is valid. To study this ques- 
tion, one has to look at systems containing smaller 
numbers of particles, since in smaller systems the 
quantum effects (in particular shell effects) which 
are neglected in the LDA, are supposed to be more 
important. 

In Fig. 2| we present the HFB (full line) and LDA 
(dashed line) results for the pairing field in the cen- 
ter of the trap, A(0), as a function of the number of 
atoms N. The calculations are done again at zero 
temperature and with a coupling constant g = — 1 
in trap units. We observe that the two calculations 
are in reasonable agreement for numbers of atoms 
greater than about 5000, which confirms the ex- 
pectation that the LDA is a valid approximation 
for systems with a large number of atoms. 

What is particularly interesting to look at in 
this figure is the region N < 3000. In this re- 
gion the HFB results clearly show the shell struc- 
ture: the pairing field becomes zero for N = 
240, 330, 440, . . . , which are the harmonic oscilla- 
tor "magic numbers". One also realizes that the 
central value of the pairing field is smaller if the 
outer shell corresponds to odd-parity states, than 



in the case where the outer shell corresponds to 
even-parity states. This can be understood eas- 
ily, since the main contribution to the pairing field 
comes from the states near the Fermi surface, and 
only s states can contribute to the pairing field at 
r — 0. Usually one expects that the LDA should 
at least reproduce the value of the pairing field if 
the fluctuations due to shell effects are averaged 
out, but our results show that the pairing field 
calculated within the LDA is systematically too 
high. This might be related to the fact that we 
are looking at the pairing field at one particular 
point (r = 0) rather than at the average gap at 
the Fermi surface, as proposed in Ref. |21j. 

When the number of atoms grows, above a value 
of about 2500 the shell structure starts to be 
washed out and gradually disappears due to the 
stronger and stronger pairing correlations. This 
happens in the region where the pairing field grows 
up to a value of about Hlo: when the pairing field 
becomes comparable with the oscillator level spac- 
ing the pairing correlations in a closed shell system 
can diffuse pairs of atoms towards the higher en- 
ergy empty shell, resulting in a non-zero pairing 
field. Globally, we observe that for N > 5000 the 
agreement between HFB and LDA is acceptable, 
even if the LDA systematically overestimates the 
value of the pairing field at the center. 

Of course, the number of particles needed for the 
validity of the LDA depends on the strength of the 
interaction; the true criterion which has to be ful- 
filled reads &lda > This criterion can even be 
applied locally, as one can see in Fig. |2 there the 
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HFB and LDA results are in perfect agreement ex- 
cept in the region of r > 5.5Z/j OJ where A becomes 
smaller than hw. 



C. Results for non-zero temperature 

Now we will discuss some results for temper- 
atures different from zero. We are particularly 
interested in the following question: Within the 
LDA, the critical temperature Tq is different at 
each point r, i.e., when the temperature increases, 
the order parameter vanishes at last in the cen- 
ter of the trap, where the local critical tempera- 
ture is the highest. In contrast to this, within the 
HFB theory, the gap and the critical temperature 
are global properties, and naively one would ex- 
pect that, as long as the temperature is below Tq, 
the pairing field extends over the whole volume of 
the system. We will see that even in cases where 
the LDA works well at zero temperature, it fails at 
non-zero temperature. On the other hand, also the 
notion that the gap vanishes globally at T — Tq, 
has to be revised in these cases. 

In Figs. El and we show the HFB and LDA 
pairing fields obtained at different temperatures, 
again for g = — 1 (in trap units) and regularization 
method (b). The chemical potentials chosen are 
/i = 32huj in Fig. [SI and [i = 40fku in Fig. corre- 
sponding to approximately 1.7 x I0 4 and 4 x 10 4 
particles, respectively. We observe that the good 
agreement obtained at zero temperature is deteri- 
orated at higher temperatures. In Fig. [5] already 
at T = 2huj/kB the LDA reproduces badly not 
only the tail of the pairing field profile, but also 
the pairing field in the central region of the trap, 
in spite of the fact that the pairing field is still 
large compared with Tilo at this temperature. The 
LDA description gets worse and worse for higher 
temperatures and results in an overestimation of 
the central pairing field and in a too drastic cut 
of the queue of the profile at large distances. Fi- 
nally, the LDA method predicts a higher critical 
temperature than the HFB one. We observed that 
Tq is equal to 3.89 (in units of Sw/fcs) for LDA 
and to 2.98 for HFB. In Fig. El the agreement is 
somewhat better. Since the critical temperature 
is higher than in the previous case, the agreement 
between LDA and HFB is maintained in a wider 
range of temperatures. Up to T — 4 one can see 
that at least the central region of the trap is well 
described by LDA. For higher temperatures, we 
observe the same kind of deterioration of the LDA 
results shown in Fig. [5] Again, the critical temper- 
ature is higher in LDA (7.08) than in HFB (5.97). 

It is evident that the LDA does not correctly de- 
scribe the phase transition in both cases. On the 
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LDA(T=2.92) 

* HFB (T=2.96) 
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FIG. 5: Pairing field A (in units of hui) as a function 
of the distance r (in units of Iho) from the center of the 
trap, for a chemical potential /j, = 32hu, corresponding 
to about 1.7 x 10 4 atoms in the trap [regularization 
method (b), cutoff Nc = 100, coupling constant g = 
— 1 in trap units]. Results obtained within numerical 
HFB calculations (symbols) are compared with LDA 
results (lines) for different temperatures T. 

15 i ■ 1 ■ 1 ■ 1 ■ 1 ■ 1 

♦ HFB (T=0) 
LDA (T=0) 
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FIG. 6: Same as Fig. O but for a chemical potential of 
fi = 40/k<J, corresponding to TV ~ 4 x 10 4 atoms in the 
trap. 



other hand, also within the HFB calculations one 
finds that with increasing temperature the pair- 
ing field becomes more and more concentrated in 
the center of the trap. Such a behavior has been 
predicted in Ref. jfl using the GL theory, the 
only assumption being that the critical tempera- 
ture is large compared with the trapping frequency, 
fcsTb ^ fiw- Let us briefly review the main re- 
sults from this theory and compare them with the 
results obtained from our HFB calculations (the 
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k F (0)\a\ 


T (0) 


T c 


5T C 


5T (GL) 


Ia 


,(GL) 


32 


0.78 


3.89 


2.98 


0.91 


1.12 


1.44 


1.23 


40 


0.91 


7.08 


5.97 


1.11 


1.29 


1.28 


0.95 



TABLE I: Comparison of results (in trap units) ob- 
tained from HFB calculations for the two cases fi = 32 
and fj, = 40 shown in Figs. and |S| [coupling con- 
stant g = 1 in trap units, regularization method (b), 
iVc = 100] and the corresponding results obtained 
from the GL theory. 



corresponding numbers are listed in Table [I]). 

In the GL theory the critical temperature Tq is 
predicted to be lower than the critical temperature 
obtained from the LDA. The difference can be 
written as 



ST C =t£ ] 



3hn /7C(3) 

; ~fczT V 48tt 2 



4fc F (0)|a| 



(34) 



where £ denotes the Riemann zeta function (C(3) = 
1.202 . . . ). In the derivation of Eq. ^ in Ref. 
the Hartree potential has been neglected. Here 
we will include the Hartree potential by using an 
effective oscillator frequency f2 > u>. Since near 
Tc the pairing field is concentrated in the center 
of the trap, we define Q by expanding the potential 
around r = 0: 

n = m^[Uo{r) + W(r))} r=0 . (35) 

Within the Thomas-Fermi approximation for the 
density profile the effective oscillator frequency can 
be written as 



2fc F (0)|q 



1 



(36) 



The estimates for STc obtained by inserting the 
numerical values for A)F(0)|a| given in Table |B into 
Eqs. (|34|l and (|36|) are very reasonable. This can be 
seen by comparing them with the STc values ob- 
tained from the HFB calculations, which are also 
listed in Table [I] If one considers that these num- 
bers can only be a rough estimate, since ksTc is 
not really very large compared with Ml, the agree- 
ment with the HFB results is very satisfying. 

Not only the critical temperature, also the shape 
of the order parameter near the critical tempera- 
ture can be obtained from the GL theory. It can 
be shown that for temperatures very close to Tc 
the pairing field has the form of a Gaussian, 



2 

A(r) = A(0)cxp(-| r 



(37) 



In contrast to the LDA result, the radius I a of this 
Gaussian is predicted to stay finite in the limit 
T — > Tc, as it is the case for the solution of the 
HFB equations. Its value is given by 



R 



mi /7C(3) 



i 



TF 



k B T\ 48tt 2 1 



4fc F (0)M 



(38) 



In Ref. |7j the quantity Rtf was defined as 
the Thomas- Fermi radius of the cloud, Rtf = 
^J2n/ (muj 2 ). Generalizing the derivation of 
Eq. (|3*8"|) to the case of a non-vanishing Hartree 
field, we see that the corresponding parameter for 
the pairing field near the center of the trap is given 
by 



Rtf 



2[fi - W(0)} 



i-MW lM0)t 



(39) 



On the other hand, the HFB pairing fields corre- 
sponding to the temperatures next to Tc shown in 
Figs.JSJandElare also perfectly fitted by Gaussians. 
As shown in Table [I] the agreement between the 
radii obtained from this fit are again in reasonable 
agreement with the radii obtained from Eqs. (1381) 
and (|39[) . The deviations are of the order of 30%, 
which is even better than one could have expected, 
since the parameter Ml/ '(ksTc) is not very small 
in the present case. 

Finally, let us look more closely at the critical 
behavior near Tc- Again, from the GL theory one 
can derive that for T — > Tc the value of the pairing 
field in the center should go to zero like 



A(0) 



1671-2^ 



T C (T C - T) . 



(40) 



As shown in Figs. and |SJ this formula is very well 
satisfied by the HFB results in both cases, fi = 32 
and fj, — 40 (in trap units) . Note that the prefactor 
in Eq. (|40|) differs from the prefactor in LDA. In 

LDA one finds for T sa T, 



c 



A 



LDA 



(0) = 



' 8tt 2 
7C(3) 



T (o )(T (o) 



T) . (41) 



The different prefactor, as well as the different crit- 
ical temperature and the finite radius of the pair- 
ing field, are due to the "kinetic" term oc AV 2 A 
in the GL energy functional, which is absent in the 
LDA and which is very important for the descrip- 
tion of the strongly r dependent pairing field near 
the critical temperature. 

As a final remark let us mention that the differ- 
ent calculations which we have compared in this 
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FIG. 7: Value of the pairing field in the center of the 
trap, A(0) (in units of ftiS), as a function of temper- 
ature T (in units of huj/ks) for a chemical potential 
fi — 32hui, corresponding to about 1.7 x 10 4 atoms in 
the trap [regularization method (b), cutoff Nc = 100, 
coupling constant g = — 1 in trap units]. Results ob- 
tained within numerical HFB calculations (symbols) 
are compared with the LDA result (dashed line) and 
with the formula 14011 obtained from the GL theory 
(solid line). 




FIG. 8: Same as Fig. [T] but for a chemical potential of 
fi = 40huj, corresponding to N ~ 4 x 10 4 atoms in the 
trap. 

paper, are all based on mean- field theory, and 
therefore do not take into account fluctuations of 
the order parameter A. It is well-known that fluc- 
tuations are very important near the phase transi- 
tion, and in particular in a situation where kp\a\ 
is not small, as it is the case here, they can lead to 
a considerable change of the critical temperature. 
Anyway, what we wanted to point out here, is that 
the LDA gives the wrong Tc as compared with a 



theory taking into account the inhomogeneity of 
the system. From this result we conclude that in 
order to have a reliable prediction of Tc for the 
trapped system, it is not sufficient to do a reliable 
calculation of Tc (even including fluctuations) for 
a homogeneous gas and then apply the LDA. 



IV. CONCLUSIONS 

In this paper we have shown a detailed compar- 
ison between HFB and LDA calculations at T = 
and at T ^ for a low density gas of superfluid 
fermionic atoms trapped by a spherical harmonic 
potential. We have used a zero-range interaction 
for the atoms and we have proposed an improve- 
ment of the regularization method adopted to re- 
move the ultraviolet divergence |8(. This improve- 
ment is a modification of aprocedure proposed for 
nuclear systems in Ref. [llj . where the Thomas- 
Fermi approximation is used in the calculation of 
the regular part of the Green's function G° reff , 
Eq. (|16fl . The use of the Thomas- Fermi approx- 
imation allows to treat systems with a large num- 
ber of atoms much easier than in the calculations 
of Ref. [8|. On the other hand, our modification 
considerably improves the convergence rate of the 
procedure with respect to the numerical cutoff. By 
using this regularization method we have observed 
that the LDA results are in quite good agreement 
with the corresponding HFB results at zero tem- 
perature and for systems with a relatively large 
number of atoms, where the shell structure effects 
are washed out. The shell effects, which are im- 
portant for small systems where the pairing field 
is smaller than the harmonic level spacing hu, can- 
not obviously be reproduced by a LDA calculation. 

For non-zero temperatures the agreement be- 
tween HFB and LDA is deteriorated even in those 
cases where it was good at T = 0. In general, LDA 
overestimates the value of the pairing field in the 
center of the trap, cuts too drastically the tail of 
the radial profile of the pairing field at large dis- 
tances, and overestimates the critical temperature 
with respect to HFB. We have verified that this 
discrepancy between the HFB and LDA results at 
T different from zero can be nicely predicted by 
using the GL theory in cases where the criti- 
cal temperature is much larger than the harmonic 
level spacing. 

In this article we considered only spherical 
traps. However, the traps used in experiments are 
usually cigar-shaped with a low longitudinal and a 
high transverse trapping frequency, lj z <C ujj_. In 
this case it is possible that the pairing field, even 
if it is larger than tuv Zl is still smaller than Huj±, 
and the LDA would probably not work. Therefore 
in principle one should also perform deformed 
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HFB calculations, but at the moment this seems 
to be numerically very difficult. On the other 
hand, as noted above, even in the case where A 
is large compared with both trapping frequencies, 
the LDA is not adequate at non-zero temperature. 
Therefore a first step to study non-spherical traps 
could be to generalize the GL theory of Ref. to 
the deformed case. 
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